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Starvation  of  polymer  electrolyte  fuel  cells  (PEFC)  takes  place,  especially  during  transients,  if  reactants 
are  consumed  in  the  fuel  cell  faster  than  they  can  be  supplied.  It  is  one  of  the  main  causes  of  aging  and 
degeneration  of  fuel  cells. 

To  prevent  oxidant  starvation  and  to  allow  for  a  dynamic  operation  of  the  fuel  cell,  the  excess  ratio 
of  oxygen  needs  to  be  adjusted  rapidly  by  increasing  the  mass  flow  into  the  cathode.  This  increase  is 
limited  by  the  inertia  of  the  actuators.  Especially  at  fast  load  changes  the  risk  of  starvation  is  high.  This 
problem  can  be  faced  by  limiting  the  dynamics  of  load  changes  or  by  decoupling  the  desired  load  from  the 
effective  load.  This  work  presents  a  decoupled  approach,  where  the  effective  fuel  cell  current  becomes  a 
controllable  variable.  Consequently,  the  control  variables  are  the  oxygen  excess  ratio,  the  fuel  cell  pressure, 
and  the  effective  current.  To  prevent  starvation  the  control  design  has  to  guarantee  that  the  oxygen  excess 
ratio  does  not  fall  beyond  a  minimum  value.  This  goal  is  achieved  by  model  predictive  control  which  is  an 
optimal  control  scheme  that  incorporates  actuator  limitations  and  state  constraints  in  the  control  design. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

An  undersupply  of  fuel  cells  leads  to  a  depletion  of  reactants 
at  the  reaction  surface  and  finally  to  fuel  cell  starvation.  Fuel  and 
oxidant  starvation  are  one  of  the  main  causes  of  short  life  and  per¬ 
formance  degradation  of  fuel  cells  [1].  When  starvation  occurs  a 
reversal  of  the  fuel  cell  voltage  can  happen  and  carbon  corrosion 
takes  place  [2,3].  Especially  during  transients,  starvation  can  occur 
if  reactants  are  consumed  in  the  fuel  cell  faster  than  they  can  be 
supplied  by  the  delivery  system.  To  prevent  starvation  corrective 
actions  as  rate  limiters  and  reference  governors  for  the  fuel  cell 
current,  as  well  as  fast  mass  flow  control  need  to  be  employed 
[4-6]. 

In  the  following  a  new  model  of  the  gas  dynamics  of  a  fuel 
cell  system,  incorporating  the  pipes  and  fittings  at  the  inlet  and 
the  outlet,  the  valve  and  the  fuel  cell  stack  is  presented.  It 
is  a  nonlinear  model  with  parameters  identified  at  a  polymer 
electrolyte  fuel  cell  system.  To  influence  the  transients  of  the 
fuel  cell  current  the  desired  load  is  decoupled  from  the  effec¬ 
tive  load.  Thereby,  the  fuel  cell  current  becomes  a  controllable 
variable. 
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To  prevent  starvation  three  actuators  -  the  valve,  the  mass 
flow  controller,  and  the  electrical  load  -  are  used  in  an  integrated 
approach.  The  presented  model  predictive  control  approach  is  a 
multivariable  control  which  incorporates  the  pressure  control,  the 
mass  flow  control,  and  the  control  of  the  fuel  cell  current  in  one 
scheme.  The  model  predictive  controller  calculates  optimal  con¬ 
trol  inputs  regarding  actuator  limitations  and  constraints  of  the 
controlled  variables. 


2.  Model  of  the  fuel  cell  system 

The  nonlinear  model  of  the  fuel  cell  system  describes  the 
dynamic  behaviour  of  the  mass  flows  and  the  pressure  drops  of 
the  gas  supply,  the  fuel  cell  stack,  the  exhaust,  and  the  outlet  valve. 
Furthermore,  a  model  of  the  mass  flow  controller  and  a  model  of 
the  electrical  load  are  incorporated. 

2.1.  Gas  dynamics 

The  model  of  the  gas  dynamics  is  based  on  the  law  of  conserva¬ 
tion  of  mass  and  the  description  of  laminar  and  turbulent  flow  [7,8]. 
The  model  does  not  account  for  spatial  dependencies.  In  analogy  to 
electronic  circuits  the  gas  dynamics  is  described  by  an  equivalent 
circuit  of  lumped  laminar  and  turbulent  flow  resistances,  Riam  and 
RtUrb»  as  well  as  mass  storing  capacitances  C  [9,10]. 
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Nomenclature 

A  system  matrix 

B  input  matrix 

C  capacity 

E  tracking  error 

/  nonlinear  function 

Hp  prediction  horizon 

J  objective  function 

k  discrete  time  index 

JCeXp  exponential  parameter  of  the  turbulent  flow  resis¬ 
tance 

Kjam  parameter  of  the  laminar  flow  resistance 
/Cturb  linear  parameter  of  the  turbulent  flow  resistance 

m  mass 

UC  vector  or  matrix  of  nonlinearities 

p  pressure 

A p  pressure  difference 

Q,  weighting  matrix 

R\am  laminar  flow  resistance 

ftturb  turbulent  flow  resistance 

R  weighting  matrix 

t  continuous  time 

T  sampling  time 

T  temperature  (with  index) 

u  input  vector 

A U  vector  of  future  changes  to  the  input  vector 

v  offset  vector 

W  mass  flow 

x  state  vector 

y  output  vector 

Y  vector  of  predicted  outputs 

Yf  free  response 

yo2  mass  fraction  of  oxygen  in  air 

Aq2  excess  ratio  of  oxygen 

Greek  letters 

a  correction  factor 

r  time  constant 

Indices 

0  standard  conditions 

c  continuous 

C  capacity 

d  discrete 

des  desired 

eff  effective 

EL  electrical  load 

FC  fuel  cell 

H20  water 

in  inlet 

max  maximum 

MFC  mass  flow  controller 

min  minimum 

02  oxygen 

out  outlet 

react  reaction 

ref  reference 

V  valve 


The  static  behaviour  of  a  component  is  determined  by  a  nonlin¬ 
ear  function 


wm 


^lam 


^turb 
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wm 
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Fig.  1.  Equivalent  circuit  of  a  component  of  the  gas  dynamics  with  laminar  and 
turbulent  flow  resistance,  R\am  and  i?tUrb>  and  capacitance  C. 


of  the  pressure  drop  Ap  over  the  component  in  dependence  on  the 
mass  flow  W  =  m  through  it,  with  the  nonlinear  flow  resistance  R 
which  is  composed  of  a  laminar  and  a  turbulent  part 

R  =  Rlam  +  Rturb  =  \  •  Warn  +  ^turb  WKFC,exp-l  ).  (2) 

The  exponential  parameter  JCexp  describes  the  nonlinear  rela¬ 
tion  of  the  pressure  drop  and  the  mass  flow  at  a  turbulent  flow 
resistance.  The  parameter  of  the  laminar  flow  resistance  J<jam  and 
the  parameters  of  the  turbulent  flow  resistance  I<turb  and  JCexp  are 
constant.  As  a  correction  factor 


To -  Pin 
^in  '  Po 


0) 


accounts  for  the  deviation  from  standard  pressure  and  standard 
temperature  at  the  inlet  of  the  component.  Therewith,  the  pressure 
difference 


A P=f(W,pin,Tin) 


(4) 


is  a  nonlinear  function  of  the  current  mass  flow  through  the  com¬ 
ponent,  the  temperature,  and  the  pressure  at  the  inlet. 

The  dynamic  behaviour  of  a  component  is  related  to  the  capa¬ 
bility  to  store  mass  in  its  volume.  With  the  differential  equation  of 
the  pressure  at  a  capacitance  C 


a tp=p= 


(5) 


and  the  equivalence  to  Kirchhoff  s  current  law  for  mass  flows 


Wm  =  WC  +  W  (6) 

with 

Ap  =  Pin-Pout  (?) 

R  R  K  J 

the  differential  equation  of  a  component 

Win  =  C  pin  +  ^  =  C  -Pin  +  Pin  RP°Ut  (8) 

results.  All  components  as  pipes,  manifolds,  fittings,  valves,  and  the 
fuel  cell  stack  itself  can  be  described  by  equation  (8)  and  the  accord¬ 
ing  equivalent  circuit  depicted  in  Fig.  1.  For  pipes,  manifolds,  and 
fittings  the  parameter  JCexp  is  equal  to  2.  For  valves,  where  the  flow 
is  purely  turbulent,  the  parameter  J<jam  =  0  disappears. 

Fig.  2  shows  the  whole  system  of  the  gas  dynamics  of  a  fuel  cell. 
All  components  in  between  the  mass  flow  controller  (MFC)  and 
the  fuel  cell  stack,  as  well  as  all  components  in  between  the  fuel 
cell  stack  and  the  outlet  valve  are  combined  to  one  component  for 
the  inlet  and  one  for  the  outlet. 
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Fig.  2.  Equivalent  circuit  of  the  fuel  cell  system  comprising  a  mass  flow  controller 
(MFC),  the  fuel  cell  stack,  inlet  and  outlet,  and  a  valve. 
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which  contributes  to  the  effective  mass  flow  in  the  fuel  cell  is, 
with  Faraday’s  law,  proportional  to  the  effective  fuel  cell  current 
JFC,  where  fcreact,o2  is  the  proportionality  constant. 

The  ratio  of  oxygen  fed  through  the  fuel  cell  WFC,o2  to  the  oxygen 
consumed  in  the  reaction  Wreact  o2 


WFc,o2 
Wreact,  > 


(20) 


Fig.  3.  Equivalent  circuit  of  the  fuel  cell  stack  with  effective  mass  flow  Weff. 


The  differential  equation  of  the  inlet  with  I<in  exp  =  2  is 

i/i/  r  n  i  Pin- PFC,in 
Win  =  On  •  Pin  +  - n - — 

*Mn 

with  the  flow  resistance 

Pin  —  '  (Pin,  lam  +  Pin,turb  ^FCJn)- 

°in 

The  differential  equation  of  the  outlet  with  I<outi  exp  —  2  is 

Wont  =  C0ut  •  PFC,out  +  ^FC’nUt - — 

Kout 

with  the  flow  resistance 


(9) 

(10) 


(11) 


Pout—  '  (Pout,  lam  +  Pout,  turb  Wout)*  (12) 

Oout 

For  the  valve  which  does  not  have  the  capability  to  store  mass 
an  algebraic  equation 

ay  •  2\pv  =<7\/  Rv  ■  Wout  =  Pv,turb  W^yfexp  (13) 


results  with  the  flow  resistance 

Rv=T.kv  bW^xp-1  (14) 

^V 

Equation  (13)  can  be  combined  with  the  differential  equation 
(11) to 

\m  r  n  i  PFC, out- Pout  Mc-s 

Wont  —  Glut  '  PFC,out  H  E -  (15) 

^Vout 

with  the  sum  of  the  flow  resistances  PVout  =  Pv  +  Pout- 

The  model  of  the  fuel  cell  stack  (Fig.  3)  differs  from  the  other 
submodels  because  of  the  electrochemical  reaction  that  is  taking 
place  and  the  mass  flows  across  the  membrane  which  are  modelled 
by  an  effective  mass  flow. 


WFC  =  VVFC  eff  +  WFc,out 

(16) 

The  differential  equation  of  the  fuel  cell  is 

m/  r  A  i  PFC,in -PFC,out 

WFc,in  =  CFC  •  PFC,in  +  p 

Kpc 

(17) 

with  the  flow  resistance 

Rpc  =  2-  .  (RFC,iam  +  RFC,turb  exp-1  ) 

orFC 

(18) 

For  the  fuel  cell  stack,  the  model  parameters  may  depend  on  the 
pressure  level,  the  humidity  and  the  current  mass  flow. 

The  structure  of  the  model  of  the  gas  dynamics  can  be  used  both, 
for  the  anode  and  the  cathode  side  of  a  fuel  cell.  The  incorporated 
parameters  were  identified  for  both  sides.  In  the  following  only  the 
cathode  side  will  be  regarded. 


2.2.  Conversion  of  oxygen 


The  fuel  cell  reaction  -  at  the  cathode  side  the  reduction  of  oxy¬ 
gen  -  couples  the  gas  dynamics  to  the  fuel  cell  current.  The  mass 
flow  of  the  reacting  oxygen 

Wreact,o2  —  kreact,02  he  (19) 


is  called  excess  ratio  and  is  defined  as  the  inverse  of  the  conversion. 
Thereby,  the  mass  flow  of  oxygen  through  the  fuel  cell  is 

WFc,o2  =  yo2  •  (WFC  -  Wfc,h2o) 

=  yo2  ■  (PFC’inRFcFC’°Ut  -  Wfc,h2o)  (21 ) 

with  the  mass  fraction  of  oxygen  in  airy02. 

For  values  ofy02  close  to  one  or  less,  not  enough  reactants  flow 
through  the  fuel  cell  to  sustain  the  reaction.  In  other  words,  starva¬ 
tion  of  fuel  cells  takes  place,  when  less  oxygen  is  fed  to  the  cathode 
as  is  consumed  in  the  reaction. 

2.3.  Models  of  the  actuators 

Besides  the  outlet  valve,  the  mass  flow  controller  and  the  elec¬ 
trical  load  are  actuators  of  the  fuel  cell  system.  The  dynamics  of  the 
mass  flow  controller  (MFC)  is  modelled  as  a  first-order  lag  element 

Wd  =  TMFC-Win  +  Win  (22) 

with  the  desired  mass  flow  Wd,  the  mass  flow  into  the  system  Win, 
and  the  time  constant  tMfc  which  describes  the  inertia  of  supply¬ 
ing  an  air  mass  flow  into  the  system.  If  models  of  compressors  or 
blowers  exist  they  can  be  used  instead. 

The  electrical  load  (EL)  is  also  modelled  as  a  first-order  lag  ele¬ 
ment 

%,d  =  rEL  •  ^FC  +  JfC  (23) 

with  the  effective  fuel  cell  current  JFC,  the  control  signal  of  the 
electrical  load  JEL  d,  and  a  much  smaller  time  constant  rEL.  By  distin¬ 
guishing  between  the  control  signal  of  the  electrical  load  rEL  d  and 
the  desired  load  current  Jd,  the  desired  and  the  effective  fuel  cell 
current  are  decoupled.  Thus,  the  current  can  be  used  as  a  control 
variable  to  avoid  starvation.  If  the  fuel  cell  is  not  connected  directly 
to  the  load,  dynamic  models  of  the  power  electronics,  e.g.  inverters, 
can  be  used  instead  of  equation  (23). 

2.4.  State  space  description 

2.4.1.  Nonlinear  continuous  state  space  description 

The  gas  dynamics  together  with  the  actuators  form  the  con¬ 
trolled  fuel  cell  system.  The  input  or  control  variables  of  the  system 

u=[wd  Rv  /FC,d]T  (24) 

are  directly  linked  to  the  mass  flow  controller,  the  outlet  valve,  and 
the  electrical  load  as  the  actuators  of  the  system.  The  state  vector 

X=  ^  Win  Pin  Ppc.in  PFC.out  ^FC  ]  (25) 

comprises  the  mass  flow  into  the  system,  the  pressure  at  the  mass 
flow  controller,  the  pressure  at  the  inlet  and  the  outlet  of  the  fuel  cell 
stack,  and  the  effective  fuel  cell  current.  The  output  or  controlled 
variables 

Y=  [V,  PFC.in  /fc]T  (26) 

are  the  excess  ratio  of  oxygen,  the  pressure  at  the  inlet  of  the  fuel 
cell  and  the  fuel  cell  current. 
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Equations  (22),  (9),  (17),  (15)  and  (23)  lead  to  a  state  space 
description 

x  =  f(x,  u)  =  Ac  x  +  Bc  u  +  A/£x(x)  +  A fCu{x)  u  (27) 

with  the  constant  system  matrix  A,  the  input  matrix  B,  the  vec¬ 
tor  of  isolated  nonlinearities  A/£x(x),  and  the  matrix  of  input  affine 
nonlinearities  AfCu(x).  The  output  equation 

y  =  A/£y(x)  (28) 

is  also  nonlinear. 


2.4.2.  Linearization 

The  vectors  A/£x(x)  and  Af£u(x)  are  linearized  by  a  first  order 
linearization  along  the  reference  trajectory  of  the  states,  element  by 
element.  For  example,  an  element  A/£x?j(xi ,  x2 )  of  the  vector  AT£u(x) 
that  depends  on  the  states  Xi  and  *2  is  approximated  by 


,  X2  )  ~  A/£x?j(x1?ref,  *2,  ref) 

dA/£X; 


+ 


dxi 


(XI  ^l,ref)  T 


dJ\fCxi 


x=xre{ 
X  (x2  -  X2  ref) 


dx2 


x=xref 


(29) 


The  matrix  A fCu{x)  is  linearized  by  a  zeroth  order  linearization 
along  the  reference  trajectory 

AfCu{x)  «  AfCu{xref).  (30) 

Therewith,  the  continuous  state  differential  equation  (27)  can 
be  rearranged  to  a  linear  state  description 

x  =  f(x,  u)  «  Acx  +  Bcu  +  vc  (31) 

with  the  system  matrix  Ac,  the  input  matrix  Bc,  and  the  offset  vector 
vc.  The  linearized  output  equation  results  to 

y  =  A fCy{x)  %Cx  +  vy  (32) 

with  the  output  matrix  C,  and  the  offset  vector  vy.  Both  offset  vectors 
merely  depend  on  the  reference  trajectory  of  the  state  vector. 


2.4.3.  Discretization 

For  the  numerical  calculation  of  the  control  variables  the  sys¬ 
tem  needs  to  be  discretized.  Thereby  it  is  assumed  that  the  input 
u  and  the  offset  vc  are  piecewise  constant  over  a  sampling  period. 
Discretization  yields  a  discrete  state  space  description 


Xfc+1  =  Ad  xk  +  Bd  uk  +  vd 


(33) 


with  the  system  matrix  Ad  =  e^7,  the  input  matrix  Bd  = 


Bc  dv  =  (Ad  -  /)  ■  A"1 
i;cdv  =  (Ad-/)-A“1 


•  Bc  and  the  offset  vector  vd 


0 

•  vc.  The  output  equation  remains  unchanged 


yk  =  Cxk+vy  (34) 

The  index  k  indicates  the  discrete  time  index  for  the  point  in 
time  tk  =  k-T,  with  the  sampling  time  T. 


3.  Model  predictive  control 

Model  predictive  control  is  a  multivariable  control  that  takes 
account  of  actuator  limitations  as  well  as  state  and  output  con¬ 
straints.  The  predictive  controller  has  an  internal  model  which  is 
used  to  predict  the  behaviour  of  the  system  depending  on  the  input 


Fig.  4.  MPC  receding  horizon,  reference  trajectory  yref,  desired  and  predicted  output 
variable,  yd  and  y. 


variables  over  a  prediction  horizon.  Thereby,  the  prediction  horizon 
Hp  is  the  number  of  time  increments  for  the  forecast. 

Figs.  4  and  5  illustrate  the  idea  of  the  receding  horizon  of  a 
single-input  and  single-output  (SISO)  system.  In  Fig.  4  the  mea¬ 
sured  output  yk  at  the  time  t  =  k-T  shows  a  difference  to  its  desired 
value  yd.  To  approach  the  desired  value  smoothly  an  exponential 
reference  trajectory 

yref  =  y<i  -  (y<i  -  y*)  •  e-  4?  (35) 

with  a  defined  time  constant  rref  is  calculated  for  the  output.  The 
task  of  the  optimal  controller  is  to  optimise  the  sequence  of  input 
changes  in  Fig.  5  so  that  the  predicted  output  y  over  a  prediction 
horizon  of  Hp  time  increments  yields  a  minimum  error  to  the  refer¬ 
ence  trajectory.  Thereby,  the  evaluation  of  the  error  can  be  limited 
to  the  time  window  te  [(Hw  +  k)T\  (Hp  +  k)T]  at  the  end  of  the  pre¬ 
diction  horizon  and  the  input  variable  can  be  left  unchanged  after 
Hu  time  increments.  If  the  optimal  trajectory  of  the  input  variable 
is  found,  the  first  input  change  is  applied  to  the  system.  For  each 
time  increment  the  prediction  horizon  is  shifted  and  the  input  is 
calculated  again. 

For  a  linear  model  the  solution  is  found  analytically.  If  the  solu¬ 
tion  does  not  fulfil  the  constraints  the  global  minimum  is  found 
solving  a  quadratic  programming  problem  by  the  active  set  method. 

3.1.  Prediction  of  the  controlled  variables 

The  system  model  in  form  of  the  linear  discrete  state  space 
description  in  equations  (33)  and  (34)  enables  a  prediction  of  the 
output  variables  by  computationally  fast  matrix- vector  operations. 
The  predicted  output  variables  over  the  prediction  horizon  are  in 
matrix- vector  form  [11,12]. 

Y=V  ■  xk+r  •  uk_!  +GX  •  vd+Gy  •  t/y+0  •  AU=Y{+&  ■  A U  (36) 

The  vector  of  the  predicted  output  variables  y  is  a  function  of  the 
current  state  vector  xfe,  the  preceding  input  vector  \xk_\ ,  the  offsets 
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vd  and  vy  and  the  future  changes  of  the  input  vector  At).  The  lengthy 
definition  of  the  matrices  4*,  T ,  Gx,  Gy  and  0  is  to  be  found  in  the 
cited  references.  The  output  is  composed  of  the  influence  of  future 
input  changes  0  •  At)  and  the  free  response  Yf  which  describes  the 
output  behaviour  in  the  case  the  input  would  not  be  changed. 

The  following  sources  of  error  lead  to  an  accumulated  error  in 
the  predicted  output: 

-  Errors  in  the  continuous  model 

-  Calculation  of  the  reference  trajectory 

-  Linearization  along  the  reference  trajectory 

-  Discretization 

-  Simplifying  assumptions  as  piecewise  constant  offsets  over  a 
sampling  period. 

For  the  fuel  cell  system  model  the  error  of  the  predicted  out¬ 
put  is  sufficiently  small  to  calculate  an  optimal  input  trajectory  as 
described  in  the  following. 

3.2.  Unconstrained  optimisation  of  the  control  variables 


The  optimum  which  is  searched  in  the  model  predictive  control 
scheme  is  characterised  by  a  minimum  error  of  the  predicted  output 
to  the  reference  trajectory,  and  additionally  a  minimum  of  changes 
of  the  input  variable  to  save  actuating  energy.  The  discrete  objective 
function  which  has  to  be  minimised  is 
Hp  Hu-1 

7—  ^  J  \9k+i  —  yref, k+i I  Iq/q  +  ^ 
i=Hw  i=0 


=  lir  — rrefllQ  +  l|A0||2  (37) 

with  the  diagonal  weighting  matrices  Qand  R.  It  can  be  expressed 
with  sums  and  in  matrix-vector  form.  As  discussed  above,  the  dif¬ 
ference  of  the  predicted  output  to  the  reference  trajectory  y  -yref 
does  not  need  to  be  penalized  at  all  times  over  the  prediction  hori¬ 
zon.  It  is  sufficient  to  penalize  only  the  error  at  the  last  Hp-Hw  time 
increments.  On  the  other  hand,  the  trajectory  of  the  input  values 
can  be  chosen  constant  after  Hu  steps.  Therewith,  less  input  moves 
need  to  be  optimised. 

Introducing  the  tracking  error 

E  =  yref  -  Yf  =  Yref  -  ( V  ■  xk+  r  ■  uk_i  +  Gx  ■  vd+ Gy  ■  vy )  (38) 

as  the  difference  between  the  reference  trajectory  and  the  free 
response  for  no  input  changes  the  objective  function  can  be  rear¬ 
ranged  [11].  To  find  the  optimal  input  changes  the  gradient  of  the 
objective  function  is  set  to  zero.  The  resulting  vector  of  the  optimal 
input  changes  can  be  calculated  analytically 

AUopt  =  0.5  •  (0TQ.  0  +  R)_1  ■  2  0tQ_E.  (39) 

To  solve  the  optimisation  numerically  it  is  formulated  as  a  least- 
squares  problem  with  the  square  roots  Sq  and  Sr  of  the  positive 
definite,  diagonal  weighting  matrices  Qand  R.  Af)opt  is  the  least- 
squares  solution  of  the  equation 


Sr 


au  = 


SqE 

0 


(40) 


which  can  be  solved  for  A0Opt  using,  e.g.  the  MATBLAB  matrix  left 
division. 


3.3.  Constrained  optimisation  of  the  control  variables 

If  the  allowed  range  of  the  input  and  output  variables  is  limited 
a  constrained  optimisation  problem  has  to  be  solved.  The  opti¬ 
misation  problem  of  the  MPC  with  inequality  constraints  has  the 


Fig.  6.  Structure  of  the  closed  loop  control  circuit. 


form  of  a  quadratic  programming  (QP)  problem.  There  are  stan¬ 
dard  algorithms  available  to  solve  QP  problems  subject  to  inequality 
constraints  as  for  example  the  active  set  method  [13,14],  which  is 
used  in  the  MATLAB-function  linprog.  If  no  constraint  are  active  the 
global  optimum  is  found  solving  the  least-squares  problem  in  equa¬ 
tion  (40).  If  one  or  more  constraints  are  active,  the  method  finds 
the  minimum  feasible  solution  of  the  QP  problem  subject  to  the 
constraints. 

4.  Control  design  for  the  fuel  cell  system 

Fig.  6  shows  the  structure  of  the  control  loop  with  the  vector  of 
the  desired  values  yd,  the  model  predictive  controller  (MPC)  which 
embraces  the  prediction  of  the  output  and  the  optimisation  of  the 
input,  and  the  nonlinear  fuel  cell  model. 

The  vector  of  the  desired  values  for  the  outputs  is 

Yd  —  [^d  PFC,d  Jd  ]  J  (41) 

where  A,d  =  2  and  pFCd  =  1.5  bar  are  constant  set  points  and  Jd  =  Jd(t) 
is  the  trajectory  of  the  desired  load  current. 

The  actuators  of  the  fuel  cell  system  are  subject  to  limitations 
which  may  be  due  to  physical  limits  or  due  to  the  balance  of  plant. 
On  the  one  hand  the  flow  resistance  of  the  valve  is  physically  lim¬ 
ited  to  a  minimum  value  determined  by  its  cross  sectional  area  and 
the  mass  flow  controller  is  only  able  to  supply  a  maximum  mass 
flow  into  the  system.  On  the  other  hand,  the  desired  mass  flow  is 
arbitrarily  limited  to  a  minimum  value  to  sustain  the  reaction  and 
to  humidify  the  membrane.  These  actuator  limitations  result  in  an 
inequality  equation  for  the  input  variables:  umin  <  uk  <  umax. 


Fig.  8.  Detail  of  the  trajectory  of  the  desired  and  the  effective  fuel  cell  current  of 
Fig.  7. 
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Fig.  9.  Trajectory  of  the  constraint  excess  ratio  A.q2  >  1. 
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Fig.  10.  Trajectory  of  the  cathode  pressure  pFC  in. 
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Figs.  7-11  show  the  response  of  the  closed  loop  to  the  load  pro¬ 
file.  As  can  be  seen  in  Fig.  7  and  in  detail  in  Fig.  8  the  MPC  reduces  the 
slope  of  the  effective  fuel  cell  current  in  comparison  to  the  desired 
current  due  to  the  exponential  reference  trajectory  and  reduces  the 
current  itself,  when  constraints  are  active  as  in  Fig.  9  the  excess  ratio. 
The  fuel  cell  current  is  not  increased  again  until  the  constraint  of  the 
excess  ratio  gets  inactive.  Subsequent,  the  effective  fuel  cell  current 
reaches  its  desired  value  without  overshooting. 

Fig.  10  shows  that  at  the  rising  fuel  cell  current  the  pressure  at 
the  fuel  cell  drops,  but  without  approaching  the  constraint.  At  the 
falling  fuel  cell  current  for  a  short  period  of  time  the  pressure  is 
higher  than  the  desired  value. 

Fig.  11  shows  the  active  actuator  limitations.  At  the  fast  load 
change  towards  higher  currents  the  desired  mass  flow  reaches  its 
maximum  value  and  keeps  it  until  the  excess  ratio  and  the  pressure 
reach  their  desired  values. 

The  MPC  successfully  performs  the  task  of  guaranteeing  a  min¬ 
imum  value  of  the  excess  ratio  and  thereby  enables  a  dynamic  fuel 
cell  operation. 

6.  Conclusion 

The  presented  model  of  the  fuel  cell  system  which  incorpo¬ 
rates  the  gas  dynamics,  the  coupling  of  the  fuel  cell  current  to 
the  mass  flow,  and  the  models  of  the  actuators  is  used  in  a  lin¬ 
earized  and  discretized  form  as  the  internal  model  of  a  predictive 
control  scheme.  The  presented  control  design  prevents  fuel  cell 
starvation  by  guaranteeing  a  minimum  value  of  the  excess  ratio 
of  oxygen  and  a  minimum  value  of  the  fuel  cell  pressure.  The  con¬ 
troller  additionally  integrates  actuator  limitations  in  its  calculation 
of  the  optimal  control  inputs.  As  the  results  show,  the  application 
of  model  predictive  control  to  the  nonlinear  model  of  the  fuel  cell 
system  increases  the  dynamics,  the  reliability  and  the  durability  of 
the  fuel  cell. 

The  presented  results  of  the  designed  controller  demonstrate 
the  applicability  of  model  predictive  control  approaches  to  non¬ 
linear  fuel  cell  models.  To  reduce  the  computational  costs  of  the 
optimisation  the  implemented  algorithm  has  to  be  transferred  into 
C-code.  If  further  reductions  are  needed  alternative  model  predic¬ 
tive  control  schemes  can  be  applied  [15,16]. 


Fig.  11.  Trajectory  of  the  mass  flow  through  the  fuel  cell  stack  WFC,  the  desired, 
limited  mass  flow  at  the  MFC  Wd,  and  the  mass  flow  into  the  system  Win. 

Since  the  target  of  the  control  design  is  to  allow  dynamic  load 
changes  without  damaging  the  fuel  cell,  the  output  variables  need 
to  be  constrained:  ymin  <  Yk  <  Ymax-  Especially  to  avoid  starvation 
the  fuel  cell  pressure  and  the  excess  ratio  of  oxygen  must  not  fall 
beyond  a  minimum  value.  Hence,  the  excess  ratio  of  oxygen  is  lim¬ 
ited  to  a  minimum  value  of  Amin  =  1.5  and  the  fuel  cell  pressure  to  a 
value  of  pmin  =  0.5  bar. 

The  control  design  and  the  simulation  of  the  nonlinear  model 
is  realised  in  MATLABISimulink  using  embedded  MATLAB  files  in 
Simulink.  Thereby,  the  sampling  time  of  the  prediction  is  set  to 
T=  1  ms  and  the  prediction  horizon  to  10  ms.  The  reference  trajec¬ 
tories  of  the  output  variables  are  generated  according  to  equation 
(35)  with  the  time  constants  =  10 ms,  rp  =  10 ms  and  t/  =  5ms  for 
the  excess  ratio,  the  pressure,  and  the  current.  As  reference  trajec¬ 
tory  for  the  state  variables  the  prediction  of  the  preceding  point  of 
time  t/<_1  are  used. 

5.  Results 

To  the  controlled  fuel  cell  system  a  rectangular  load  profile  as 
depicted  in  Fig.  7  is  applied. 
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